Approved  for  public  release;  distribution  is  unlimited. 


Monte  Carlo  Simulation  of  Seismic  Location  Errors 

for  Moving  Vehicles 

October  4,  2001 

Roy  J.  Greenfield+  and  M.  L.  Moran* 

+Penn  State  University,  State  College,  PA 
*USACE  ERDC-CRREL 


ABSTRACT 


A  method  was  developed  to  predict  the  distribution  of  seismic  source  mislocation  errors  for  tracked 
vehicles  using  bearing  and  range  (r)  location  statistics  obtained  from  field  data.  The  method  is  of  use  for 
1)  guiding  the  design  and  deployment  of  seismic  sensor  networks  and  2)  statistically  modeling  seismic 
network  tracking  performance. 

In  the  present  work,  a  seismic  array  is  composed  of  one  or  multiple  seismic  subarrays  each  of  which 
makes  estimates  of  bearing  and  range.  The  heal  ing  and  range  estimates  are  taken  to  be  normally 
distributed  around  the  true  target  location  with  errors  cb  and  or  respectively.  When  these  0’s  are  used  in 
our  Monte  Carlo  simulation,  we  chose  a  true  source  point,  and  network  geometry  (array  locations).  Using 
these  source  and  sensor  positions,  a  set  of  true  bearing  and  r  are  computed.  The  Monte  Carlo  simulation 
is  done  with  100  trials,  in  each  trial  a  random  error  is  added  to  the  true  bearing  and  r  for  each  array 
according  to  the  observed  0.  For  each  of  the  100  bearing  and  r  Monte  Carlo  “measurements”,  a  weighted 
least  squares  solution  for  the  East  and  North  value  of  the  source  is  found.  A  statistical  fit  for  all  100  trials 
is  used  to  find  the  90%  confidence  location  ellipse  surrounding  the  true  source  point. 

Two  types  of  displays  are  presented.  The  first  fixes  0b  and  0r  and  makes  a  map  that  displays  the  error 
ellipse  axes  and  orientation  as  a  function  of  map  position.  The  second  uses  a  fixed  map  position  and 
displays  the  error  ellipse  axes  and  orientation  as  a  function  of  0b  and  0r.  This  allows  convenient 
analysis  of  how  the  location  error  changes  with  changes  of  the  array  measurement  accuracy.  The  later 
format  is  particularly  suited  for  evaluating  the  impact  of  seismic  senor  configuration  on  network  tracking 
performance. 

Results  are  presented  for  a  number  of  different  network  configurations  representative  of  battlefield 
situations.  Both  single  array  and  multiple  array  networks  are  considered.  The  values  used  for  0b  and  0r 
in  these  simulations  were  chosen  by  examining  results  from  a  number  of  past  field  test  data  sets. 

Using  field  data  it  was  found  that  subarrays  with  1 1  or  12  vertical  seismic  sensors  could  get  bearing  with 
RMS  errors  on  the  order  of  5°.  Subarrays  with  only  3  or  4  sensors  and  diameters  of  5  to  10  m  gave 
similar  errors.  For  small  subarrays  with  3  sensors  in  an  equilateral  triangle  with  legs  of  1  m,  scatter  was 
below  10°  but  a  significant  bias  existed;  the  bias  might  be  removed  by  a  calibration  procedure. 

Range  was  found  from  amplitude  measurements  at  five  sites.  A  model  appropriate  to  the  amplitude  decay 
of  a  surface  wave  on  lossy  medium  was  fit  to  data.  Several  band  passes  were  examined.  RMS  Range 
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errors  for  the  five  sites  varied  between  45  and  145  m.  A  curve  was  produced  comparing  amplitude 
versus  range  for  the  five  sites.  This  demonstrates  that  a  site’s  amplitude  versus  distance  must  be 
calibrated  to  estimate  distance  from  amplitude. 


1.  INTRODUCTION 


The  major  effort  has  been  towards  determining  how  accurately  a  continuous  source  of  seismic 
signal,  such  as  a  tank,  can  be  located  from  it’s  seismic  emissions.  The  location  of  a  source  is  a  key  piece 
of  information  for  monitoring  hostile  movements. 

A  seismic  location  system  is  considered  that  is  composed  of  a  group  of  subarrays  Each  subarray  has  a 
number  of  individual  seismic  sensors.  A  typical  subarray  will  have  a  diameter  of  a  few  meters  to  40 
meters.  .  All  the  subarrays  will  be  termed  an  array. 

This  paper  is  divide  into  3  sections:  1)  RANGE  ESTIMATION,  2)  BEARING  ESTIMATION,  and  3) 
LOCATION  ACCURACY.  In  the  first  two  sections  estimates  are  made  of  the  range  and  bearing 
accuracy  that  can  be  obtained  by  processing  past  field  test  data.  In  the  third  section  a  method  is 
developed  an  applied  to  estimate  source  location  errors. 

2.  RANGE  ESTIMATION 


2.1.  Method 

The  amplitude  of  a  seismic  signal  decays  because  of  geometric  spreading  and  attenuation.  The 
attenuation  is  due  to  inelastic  losses  and  to  scattering  due  to  inhomogenaity.  For  a  Rayleigh  surface  wave 
the  amplitude  dependence  on  distance  is  of  the  form 


A  =  [reference  amplitude] [geometric  spreading] [attenuation] 


where 


A0  =  reference  amplitude 


r  =  range 
Q  =  quality  factor 
c  =  wave  velocity 
/  =  frequency 


A  linear  least  squares  fit,  for  a  and  b,  is  done  in  the  linear  form 


PdB  +101og(r)  =  a  +  br 
where 

a  =  101og(i4^ ) 


b  =  10  log 


e 


Qc 


=  -20 7Uf  log  (<?)  /  Qc 


V  J 

P  =  power  and  PdB  is  the  power  in  dB. 

^,=101og(P) 


The  values  of  a  and  b  are  found  by  a  linear  least  squares  fit  to  minimize  the  form: 


E2 (a,b) 
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1=1 

Numb.Obs.  ^ 

X  [Ct+lOlog  (r)-(a  +  br)] 


2.2.  Results 


Data  from  five  of  field  tests  were  fit  and  used  to  estimate  error  in  range  estimation;  the  data  designations 
are  given  in  Table  1 .  An  example  of  a  range  fit  is  shown  in  Figure  1 .  Fits  were  made  over  several 
different  band  passes. 


Table  1 .  Data  used  in  range  fits  and  in  bearing  fits. 


Designation 

Data  Used 

Yuma.  AZ 

Yuma  Proving  Ground,  Smart  Weapons  Test  Range  9/14/200 

Four  wheel  Drive;  File  14, 

Aberdeen,  MD, 
Site  1 

June  11,  1996,  10c  runs,  Piston  Tank  ; 

10:18  34:42  53:55  58:64  72:79  92:97  105:118  %  file  vector 

Ft.  Greely,  AK, 
Site  1 

1/27/1997,  ,  Piston  Tank;  34:42  53:64  %  file  vector 

Ft.  Greely,  AK, 
Site  2 

Dec  11,  1997  ;  File  56,  ,  Piston  Tank 

Aberdeen,  MD, 

Site  2 

10/28/97  File84;  File  56 

Abdem96A4  ;  Pow./Hz 


Range  (m) 


Figure  1.  Fit  to  amplitude(Power/Hz)  versus  range  data  for  Aberdeen,  MD,  Site  1.  Band  pass, 
slope,  and  Standard  Deviation  (std)  are  given  on  each  panel. 


After  the  a  and  b  coefficients  are  found  the  range  can  be  estimated  from  a  power  level  measurement. 
Then  for  a  group  of  observed  power  levels  and  GPS  based  range  measurements  (assumed  to  be  correct) 
the  standard  deviation  of  the  error  in  the  estimated  range  is  found.  The  error  in  range  is  defined  as  the 
predicted  range  minus  the  GPS  range.  The  error  in  range  that  corresponds  to  the  fit  shown  in  Figure  1  is 
shown  in  Figure  1 .  Standard  deviations  of  the  range  fit  estimates  for  5  sites  are  shown  in  Table  2. 


Abdern96A4  ;  Pow./Hz 


Range  (m) 


Figure  2.  Range  error  versus  range  for  data  for  Aberdeen,  MD,  Site  1  (see  Table  1).  Band  pass 
and  RMS  error  are  given  on  each  panel. 


Table  2.  Compilation  of  results  for  errors  in  range  estimation.  10  to  20  Hz  band 


Data  Set  number 

Location 

Standard 
Deviation  (m) 

Range  of  Fit 
(m) 

1 

Ft.  Greely  Site  1 

145 

250-1500 

2 

Ft  Greely  Site  2 

37 

200-800 

3 

Aberdeen  Stitel 

80 

350-1100 

4 

Aberdeen  Site  2 

64 

110-350 

5 

Yuma.  AZ. 

40 

100-320 

2.3.  Comparison  of  Amplitude  for  Different  Sites 

The  fit  to  data  for  5  different  sites  was  complied  on  one  plot,  shown  in  Figure  2,  to  show  the  strong 
influence  that  geologic  structure  has  on  the  amplitude  versus  distance  curve.  For  each  site  the  curve 
represents  the  fit  to  the  data.  The  range  of  the  plot  indicates  the  range  of  the  data  used  in  the  fit.  There 
are  large  differences  in  amplitude  levels  and  in  the  attenuation  rate  for  the  different  sites.  This  indicates 
that  some  form  of  calibration  must  be  done  at  a  site  if  amplitude  measurements  are  to  be  used  to  estimate 
range.  The  coefficients  a  and  b  are  given  in  Tables  3  and  4. 
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3  Amplitude  versus  distance  for  5  sites.  Shows  strong  effect  of  geology  on  amplitude.  Lines 
extend  over  range  of  data  used  in  fit.  All  sources  were  piston  tanks  except  Yuma  which  was  a 
Sports  Utility  Vehicle. 
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Table  3.  Coefficient  b. 


Wide  Band 

5-20  Hz 

10-20  Hz 

20-40  Hz 

Yuma 

-0.0787 

-0.0717 

-0.0732 

-0.0926 

Aberdeen,  Site  1 

-0.016282 

-0.015713 

-0.018326 

-0.023465 

Ft.  Greely,  Site  1 

-0.022511 

-0.024814 

-0.028239 

-0.021444 

Ft.  Greely,  Site  2 

0.0913 

-0.0669 

-0.0667 

-0.0978 

Aberdeen,  Site  2 

-0.044983 

-0.033857 

-0.048274 

-0.089925 

Table  4.  Coefficient  a. 


Wide  Band 

5-20  Hz 

10-20  Hz 

20-40  Hz 

Yuma 

-120.8522 

-119.1978 

-118.305 

-113.7211 

Aberdeen,  Site  1 

-93.9392 

-81.7267 

-82.4574 

-96.7581 

Ft.  Greely,  Site  1 

-134.7292 

-124.7725 

-122.3322 

-121.1379 

Ft.  Greely,  Site  2 

-63.0317 

-68.164 

-67.1417 

-48.3767 

Aberdeen,  Site  2 

-95.5133 

-87.9894 

-86.1321 

-83.2409 

2.4.  Conclusions  on  Range  Accuracy 

•  Data  from  5  sites  was  analyzed. 

•  A  form  suited  for  surface  wave  amplitude  versus  distance  with  attenuation  was  used  in  the  fits. 

•  A  plot  comparing  the  amplitude  versus  range  for  the  five  sites  was  given. 

•  Range  error  RMS  for  the  sites  varied  between  40  and  145  m. 


3.  BEARING  ESTIMATION 

A  study  of  bearing  results  was  done  to  determine  what  accuracy  can  be  expected  from  well  populated 
seismic  subarrays.  This  gives  guidance  on  the  variance  to  use  in  the  simulation  results  section.  In  past 
work  (e.g.  Greenfield  and  Moran,  1997)  looking  at  bearing  determination  for  moving  vehicles  the 
subarrays  were  typically  12  to  30  m  in  diameter  and  had  12  geophones.  A  second  purpose  of  this  study 
is  to  look  at  the  use  of  small  arrays  with  limited  numbers  of  sensors  to  see  to  what  extent  bearing 
determination  deteriorates  these  smaller  subarrays  are  used.  To  do  this  processing  is  done  with  3  or  4 
sensors  selected  from  the  full  subarray.  The  past  subarray  collections  were  not  done  with  this  application 
in  mind  so  the  best  small  subarrays  that  could  be  formed  from  available  data  were  used. 


3.1.  Method 


Normal  Frequency  Wavenumber  (fk)  (e.g.  Lacoss  et.  al.,  1969)  processing  has  given  the  most  reliable  and 
consistent  results  for  determining  bearing  with  subarrays  and  is  used  in  all  results  presented. 


3.2.  Results 


An  example  of  bearing  errors  made  with  the  three  subarrays  A,  B,  and  C,  shown  in  Figure  4,  is  shown  in 
Figure  5.  Subarray  A  was  a  large  aperture,  40  m,  subarray  with  12  vertical  geophones.  It  has  been 
found  in  the  past  that  this  subarray  gives  accurate  results.  The  smaller  two  subarrays  had  diameters  up  on 
the  order  of  15  to  20  m.  The  range  to  target  was  up  to  approximately  1,400  m.  The  bearing  determination 
was  done  at  18  Hz.  The  mean  phase  velocity,  estimated  with  subarray  A,  was  1850  m/s.  A  plot  of  the 
bearings  determined  by  the  subarrays  are  given  in  Figure  3;  Figure  4  gives  the  bearing  errors.  Table  5 
gives  the  errors  found  for  the  three  subarrays.  Errors  >12°  were  excluded.  The  bearing  accuracy  of  the 
three  subarrays  was  similar,  with  RMS  errors  of  under  5°.  The  largest  errors  occur  at  the  longest  ranges 
which  may  be  a  signal  to  noise  problem.  Except  for  the  bearings  for  ranges  near  1500  m,  there  is  no 
clear  increase  of  error  with  distance.  Bearing  estimates  were  also  done  for  data  From  Aberdeen  site  1. 

The  results  are  presented  in  Table  6. 


Sub.  A  Sub.  B  Sub.  C 


x  (m)  x  (m)  x  (m) 


Figure  4.  Subarrays  for  Ft.  Greely  data.  A)  is  a  typical  subarray.  B)  and  C)  show  small  arrays 
formed  by  selected  elements  of  A) 


File  Number 

Figure  5.  Top  Ft.  Greely,  Jan  27, 1997  bearings  determined  with  3  subarrays.  See  Figure  4  to 
identify  symbols.  Small  black  dots  indicate  GPs  bearing.  Bottom  shows  range  to  tank  source. 


Table  5.  Errors  for  Ft.  Greely  data  of  Figure  6. 


subarray 

mean  error  (°) 

Standard  Deviation  (°) 

RMS  Error 

(°) 

A 

1.0934 

4.6006 

4.6855 

B 

0.7559 

4.8745 

4.8812 

C 

1.1633 

4.4129 

4.517 

Table  6.  Errors  for  Aberdeen  site  1. 


subarray 

Diameter 

(m) 

Number  of 

geophones 

mean 

(°) 

Std 

(°) 

RMS 

(°) 

A 

20 

12 

-1.5468 

5.3389 

5.5349 

B 

10 

4 

-2.897 

6.8006 

7.3631 

C 

10 

4 

-3.1525 

6.9962 

7.6443 

D 

5 

3 

6.2212 

8.839 

10.7744 

3.3.  Conclusions  on  Bearing  estimation 


•  Data  from  2  sites  was  used  to  estimate  bearing  error. 

•  Bearings  obtained  were  compared  to  GPS. 

•  Subarrays  with  approximately  12  sensors  gave  RMS  bearing  error  on  the  order  of  5°. 

•  Small  subarrays  with  3  and  4  sensors  with  diameters  of  5  to  10  m  gave  similar  errors. 

•  Data  from  a  third  site,  Yuma  AZ  was  used.  Very  small  subarrays  with  3  sensors  and  diameters  on 
the  order  of  1  m  gave  constant  bearings,  but  in  some  cases  the  bearings  had  large  biases.  These 
biases  might  be  correctable  by  careful  planting,  or  sensor  calibration. 


4.  LOCATION  ACCURACY 


4.1.  Method 

The  location  estimate  (x,  y  coordinates)  is  found  from  the  estimated  bearing  and  range  from  a  set  of 
subarrays.  The  criteria  used  is  the  non-linear  minimization  of  the  weighted  least  squares  error.  The 
method  is  called  Chi-Square  Fitting  (Press  et  al,  1992,  p  653-655).  The  quantity  minimized  is 


X3(^.3')  =  I 


'  P~P(x,y)^ 

h 

2 

N  . 

+II 

f  r-r(x,y)^ 

l  ) 

i=l 

l  <*i  J 

where 

x,  y  =  location  positions 

X2  (x,  y)  =  CHI  squared  function  to  be  minimized 
N  =  number  of  stations 

j3  =  observed  bearing  ( 1 ) 

j3(x,y)  =  predicted  bearing 

( 7 f  =  square  root  of  the  variance  of  bearing  i  measurement 
r  =  observed  range 
r(x,  y)  =  predicted  range 

a'  =  square  root  of  the  variance  of  range  i  measurement 


The  minimization  is  implemented  to  find  x  and  y  from  a  set  of  observations  using  the  Matlab  function 

Isqnonlin.m 


4.2.  Simulation  Procedure  for  a  error  confidence  ellipse  map 

An  checkerboard  grid  of  true  source  points  is  chosen  A  location  determination  error  ellipse  is  simulated 
with  the  steps  described  in  Table  7.  Then  the  error  ellipse  is  plotted  at  each  true  source  point.  A  map  is 
made  for  the  grid  of  source  points.  The  program  MapE.m  does  this  procedure. 


Table  5.  Procedure  for  producing  an  error  ellipse 


1 

select  array  configuration 

2 

select  true  x  and  y 

3 

select  a  for  range  and  bearing 

4 

generate  true  bearing  and  range 

5 

use  a  Normal  Distribution  random  number  generator  for  error  in 
bearing 

6 

use  a  Normal  Distribution  random  number  generator  for  error  in 
range 

7 

add  error  in  bearing  to  true  bearing  to  get  observed  bearing 

8 

add  error  in  range  to  true  bearing  to  get  observed  range 

9 

fit  for  x  and  y  using  X2  criteria  ( the  variances  used  in  the  random 
number  generation  are  used  in  the  fit) 

10 

repeat  steps  5  to  9  M  times 

11 

fit  the  M  locations  found  with  an  error  confidence  ellipse 

12 

plot  the  axes  of  the  error  confidence  ellipse  on  map 

Figure  6  illustrates  how  the  scatter  of  estimates  of  the  source  point  are  fit  by  an  ellipse  and  the  axes  of  the 
ellipse  are  plotted  for  each  true  source  point.  A  color  code  is  used  to  allow  a  wide  range  of  ellipse  sizes  to 
be  plotted  on  the  same  figure.  In  black  and  white  the  color  code  does  not  show  up  so  a  few  ellipse  axes 
have  been  labeled. 


Figure  6.  Formation  of  error  ellipse  maps.  Green  dots  in  left  panel  show  individual  locations 
scattered  around  true  location.  Small  figure  shows  fit  of  ellipse  to  the  individual  locations.  Semi- 


major  and  semi-minor  axis  of  ellipse  is  drawn  at  true  location.  Thin  black  line  is  length  of  the 
semi-major  axis  of  maximum  for  the  color. 


For  MapE.m  the  values  of  01  and  0r  can  be  functions  of  range.  In  the  present  implementation  <7b 
is  treated  as  constant  with  range.  <j'  is  taken  in  the  form 

r=\  ^ 0  r  <  r0 

G  ~  {°o  [l+(r  -  ro)‘^,/]  r>=ro 

A  program  was  also  written  to  explore  how  location  accuracy  at  a  fixed  map  point  varied  with  errors.  For 
this  program  VaryErr.m  the  values  of  0  and  0  are  given  as  single  values  for  the  map  point. 


4.3.  Results 

An  example  of  an  error  map  is  given  in  Figure  7.  On  these  maps  color  is  used  as  a  scaling  key,  so  in 
black  and  white  the  length  of  some  of  the  axes  are  not  clear,  thus  a  few  numbers  have  been  added. 


nw  1(B  *  W  HO  IW  P  M  ffl 


Figure  7.  Error  ellipse  map  for  three  subarrays  (circles) .  r0=80  m;  Co=  40  m;  S=.01  (1/m);  OBear. 

=  8°.  Circles  show  the  5  subarray  locations.  The  numbers  indicates  semi-major  error  ellipse  axis. 
Horizontal  axis  is  East,  Vertical  North. 


Figure  8  show  the  way  the  error  ellipse  varies  with 
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Figure  8.  Error  ellipse  dependence  on  (Rearing  and  aRange  •  Same  three  subarray  configuration  as 
previous  figure.  Source  location  at  200  East,  300  North.  Horizontal  axis  is  East,  Vertical  North. 
The  numbers  indicates  semi-major  error  ellipse  axis. 


4.4.  Conclusions  on  statistical  modeling  of  source  location  error 


•  A  weighted  least  squares  method  was  developed  to  fit  for  source  location  using  bearing  and  range 
estimates  from  one  or  more  subarrays. 

•  The  method  was  applied  to  synthetic  data,  with  known  error  distributions. 

•  Maps  giving  90%  error  ellipses  were  produced  for  a  number  of  reasonable  array  configurations. 
Errors  variances  used  were  consistent  with  values  found  from  past  field  test  data. 

•  A  presentation  was  developed  to  illustrate  the  dependence  of  the  90%  error  ellipses  on  the  range 
and  bearing  error  variances. 
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